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Abstract 

We present molecular dynamics simulations of solid NaN02 using pair potentials with the rigid- 
ion model. The crystal potential surface is calculated by using an a priori method which integrates 
the ah initio calculations with the Gordon-Kim electron gas theory. This approach is carefully 
examined by using different population analysis methods and comparing the intermolecular inter- 
actions resulting from this approach with those from the ah initio Hartree-Fock calculations. Our 
numerics shows that the ferroelectric-paraelectric phase transition in solid NaN02 is triggered by 
rotation of the nitrite ions around the crystallographical c axis, in agreement with recent X-ray 
experiments [Gohda et al, Phys. Rev. B 63, 14101 (2000)]. The crystal-field effects on the nitrite 
ion are also addressed. Remarkable internal charge-transfer effect is found. 

PACS numbers: 64.60. Cn,61.43.Bn,64.70.Pf 
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FIG. 1: Crystal structure of NaN02 in the ferroelectric phase. 
I. INTRODUCTION 

Sodium nitrite is a ferroelectric at room temperature. It has the orthorhombic structure, 
space group C|° — /m2m, with the dipole vector of the V-shaped nitrite anions aligned paral- 
lel to the crystallographic b direction, as shown in Fig.^ The ferroelectric-paraelectric phase 
transition takes place at about 437 K, where the high temperature phase is orthorhombic, 
space group Dlf^ — Immm, with the dipoles disordered with respect to the b axis. In a narrow 
temperature range from 435.5 K to 437 K, there exists an incommensurate antiferroelectric 
phase. The melting temperature is 550 K. Distinguished from displacive ferroelectrics in 
which the ferroelectric transition is driven by soft phonon modes, NaN02 offers a model 
system for research of the order- disorder structural phase transition and any associated 
ferroelectric instability.^*^ 

Extensive experimental work on NaN02 has been devoted to probing the mechanism of 
the NO^ polarization reversal that triggers the order- disorder transition. The majority of 
studies support the c-axis rotation model, but there were also results favoring the a-axis ro- 
tation model.— Recently, refined X-ray studies over a wide temperature range reinforced the 
c-axis rotation model.-^ On the theoretical side, the microscopic model calculations done by 
Ehrhardt and Michel supported the c-axis rotation mechanism,^ whereas mixed double rota- 



2 



tions around the a-axis and the c-axis was suggested by Kinase and Takahashi.— It has long 
been desirable to apply computer molecular dynamics (MD) simulations to NaN02 in order 
to achieve unambiguous understanding of the polarization reversal mechanism. Earlier MD 
simulations with empirical Born-Mayer pair potentials detected the c-axis rotation in above- 
room-temperature NaN02-^ ^''' ' ' - ° Unfortunately, the low-temperature structure produced by 
those simulations was antiferroelectric and apparently disagreed with the experimental ob- 
servations. 

Lu and Hardy pointed out that the overall phase behavior of NaN02 could be simu- 
lated by using an a priori approach to construct the crystal potential surface (PES).— The 
Lu-Hardy (LH) approach was originally designed to deal with molecular crystals such as 
K2Se04, where exists a mix of bonding types, that is, the intermolecular interactions are 
mostly ionic, but the constituent atoms in a molecule (Se04~ in K2Se04) bond covalently. 
In the LH approach, the intra-molecule interactions were treated by applying the ah ini- 
tio self-consistent field method to the gas-phase molecules, while the intermolecular pair 
potentials were computed within the Gordon-Kim (GK) electron gas theory.— The crux of 
their application of the GK theory is how to partition the ah initio molecular charge density 
between the constituent atoms. Since there is no unique way to separate the charge density 
of a highly covalently bonded molecule, Lu and Hardy suggested equal separation in a spirit 
similar to the MuUiken population analysis (MPA). By using this atomic-level method, we 
could successfully describe the phase transitions in fiuoroperovskites,— and ionic crystals 
with polyatomic molecules including SeO^,i^ ClOj,^ S0^,^ ^iOt ^ and NO3 
Note that the MPA happens to preserve the (zero) dipole moment of these molecules. 

However, several problems appear when we moved on to deal with NaN02 where the 
NO^ radical has nonzero dipole moment and stronger chemical bonding. First, it is well 
known that the MPA, while certainly the most widely employed, is also somewhat arbitrary 
and the most criticized.— In particular, the MPA overestimates the dipole moment of the 
free NO^ ion by about 120%. Other difficulties involved the free- ion approximation. Unlike 
in monatomic ionic crystals, there may exist considerable internal charge-transfer effects 
in molecular ionic crystals. Electronic band structure calculations^^ indicated that within 
a nitrite entity, the nitrogen atom and two oxygen atoms bond covalently, leading to high 
charge transferability between these constituent atoms. Therefore, in solid NaN02 the NO2 
group will feel different crystal-field environments as it rotates and responds by redistributing 



3 



the charge density among its three constituent atoms. 

Our goals in this paper are twofold. First, we show that our atomistic level simulation 
methods involving pair potentials with the rigid-ion model is capable of correctly describing 
the phase behavior of NaN02. Second, we systematically examine the LH approach to 
understand the reason why it works so well in molecular ionic crystal systems by the following 
steps: (i) we develop another population analysis method that preserves the molecular dipole 
moment by directly fitting the ab initio charge density of a molecule; (ii) we carry out 
ab initio Hartree-Fock (HF) calculations of the intermolecular interactions and find that 
the pair potentials from the rigid-ion model can correctly reproduce the ab initio results; 
(iii) we investigate the crystal-field effects on the NOg ion by embedding the ion (and 
its first shell of neighbors) in a lattice of point charges and find a remarkable internal 
charge-transfer effect.— Several MD simulations based on these modifications of the LH 
approach are also performed. The ferroelectric-paraelectric transition triggered by the c- 
axis rotation of the nitrite ions is observed in all versions of the LH approach. However, 
the transition temperatures predicted by these simulations are lower than the experimental 
values. Furthermore, the transition temperatures obtained from the original version are 
higher than those predicted by modified versions and closer to the experimental values. 
After careful examination, we notice that in the LH approach, the NO^ dipole moments 
were generally overestimated by about 120%, which reinforces the ground state. This implies 
that the crystal structure of NaN02 is stabilized by the anion polarization effects. Thus, we 
conclude polarization effects are particularly important for the quantitative study of NaN02. 

This paper is organized as follows. Section II describes the method we used to obtain 
the intermolecular interactions. Section III analyzes the resulting intermolecular potentials 
in comparison with those obtained from ab initio calculations. Section IV presents the 
results of our MD simulations. The crystal-field effects on NO2 are discussed in Section V. 
Concluding remarks are made in Section VI. 

II. INTERMOLECULAR INTERACTIONS 

Our MD simulation technique method originates from the GK model for simple ionic 
crystals such as alkali halides, assuming that (molecular) ions in a crystal environment may 
be described as free ions .^^'^^'^^ Then it was extended to deal with molecular ionic crystals 



4 



like K2Se04, in which strong intramolecular covalency exists . The main idea is that the 
molecular ion (SeOj in K2Se04) is treated as a single entity, and intramolecular and inter- 
molecular interactions are treated separately: first we perform ab initio quantum chemistry 
calculations for the whole molecular ion, to obtain the optimized structure, the force con- 
stants, and the whole electron density p(r). The intramolecular interactions are described 
by force constants within the harmonic approximation. As for the intermolecular interac- 
tions, we have to carry out electron population analysis to separate p(r) onto each individual 
atom in the molecular ion, then use the Gordon-Kim electron gas model to calculate the 
intermolecular pair potentials. This approach provides a parameter-free description for the 
crystal potential-energy surfaces, which allow structural relaxation, MD simulation, and 
lattice dynamics calculations. 

In calculating the intermolecular forces, there are three major approximations as discussed 
in the following: 

(1) We assume that the geometries and electronic densities of the separate ions remain 
unchanged once they have been obtained under given circumstance, such as in the equilib- 
rium state of the gas or crystal phases. This approximation is the fundamental basis for 
the GK electron gas theory. General speaking, we found that in an ionic crystal there is no 
strong chemical bond between ions, hence this approximation is reasonable. 

(2) When dealing with the intermolecular interaction, we assume that the charge density 
of a rigid ion can be separated into its atomic constituents. 

(3) We assume that the crystal potential energy is composed of the intermolecular and 
intramolecular interaction, where the intramolecular interaction is expressed in terms of 
force fields and the intermolecular interaction is a sum of interatomic pair potentials. 

Atomistic level simulations utilizing pair potentials and the rigid-ion model have great 
success in describing many ionic systems.— We showed that this scheme can correctly de- 
scribe the phase transition behaviors of alkali halide fiuoperovskites,-^ and molecular crystals 
with tetrahedrali^iiiJ^iiSiil and equilateral triangula r ^^i^^'^'-' radicals. However, for NaN02 in 
which NO^ has only a two fold symmetrical axis, the results were less satisfactory.— Note 
that the mean HF polarizability of NO2 , 14.156 (atomic units), calculated with the D95* 
basis,«2fi. is much higher than that of Na"*", 0.343 (atomic units), calculated with the 6-31* 
basis. Therefore, in solid NaN02, the rotation of NO2 in the crystal field will induce charge 
redistribution within the molecule. Hence this dynamic effect may invalidate the rigid ion 
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approximation. In this paper we shall perform HF calculations for various geometries to 
verify this scenario. 



A. Pairwise additive approximation 

In the GK model, we evaluate the interaction between two molecules based on the elec- 
tron density,^' which is approximated as the sum of component densities taken from HF 
calculations. That is, if pA and pb are the component densities, then the total density 
is Pab — Pa^ Pb- Whereas, interaction potential is computed as the sum of four terms: 
Coulombic, kinetic, exchange, and correlation energies which are expressed in terms of the 
charge densities. 

Therefore, suppose the A and B molecules are made up of M and atoms, respectively, 
then the Coulombic interaction between them is 



J J |ri — -^-i J |r2 — rt,A,i 



i=l 

E^bJ'^'.U^^ + EEtS^^^. (1) 



where Za^i^ ^Bj-, R-a,?. and R_B,j are the nuclear charges and coordinations of the z-th atom 
in the A molecule and j-th atom in the B molecule, respectively. This potential energy can 
be split into two parts: first the long-range part, 

^c = flJ2 ^ ^ ^"^^^^^^^^^^^^'^ -IpBi^2)dr2] 

i=lj=l \^A,i -'RbjI 

and the short-range part 

= Vc- V^. (3) 

Eq. Q is essentially the electrostatic interaction energy when the charge densities of the 
molecules are distributed as point charges on the constituent atoms, which is known as the 
Madelung potential energy. 

The non-Coulombic energy terms are expressed in the uniform electron gas formula, 

Vi = J dr[p^J^{r)e^{pAB) - PA(^)e^(PA) - PBi^hipB)]: (4) 
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where ei(p) is one of the energy functional for the kinetic, exchange, and correlation 
interactions.— Note that Eq. ^ is not composed of pair potentials. In order to obtain 
the effective pairwise potentials, we approximate Eq. (jH) using 



where p^,„ = Pm + Pn- Pm ^'^d p„ are the charge densities of individual atoms in the A and 
B molecules, respectively, which are obtained by a population analysis as described in the 
next subsection. 

Even though the non-Coulombic forces as determined by Eq. are not strictly additive, 
the above approximation appears to be adequate except at very short distances. As pointed 
out by Waldman and Gordon,— the main reason as to why this approximation is valid is 
because the Coulombic force, the largest contribution to the potentials, is additive. Based 
on our calculations, we find additivity of Vi holds only to within about 50%; however, the 
overlap contribution to the electrostatic energy dominates Vi and renders additivity to within 
10%. One final remark is in order, for the sake of simplifying the two-electron integral in 
Eq. ((T)), the charge densities, p„ and p„, are taken as its spherical average. As a result, the 
Coulombic interaction is not exactly evaluated. Nevertheless, as we shall show in Figs. El 
and 01 this error is compensated by those due to the pairwise additive approximation. 

To summarize this subsection, we demonstrated that it is possible to analytically express 
the intermolecular potentials + Vq + Vi using Eqs. 0, © and (0) once the charge 
density of each individual atom is obtained by an electronic population analysis. In the next 
subsection, we shall present further analysis on the charge density. 

B. Electronic Population Analysis 

In this subsection, we discuss the ways to separate the electron density p(r) of a molecule 
into its atomic constituents. Suppose the molecule consists of M atoms, then the wave func- 
tion of the molecule ip{r) can be written as a linear superposition of atomic wave functions 
(p{r — Rj), 2 = 1,2,..., M, centered at each atom. 



^^-Y.Y. c?r[p„„(r)e,(p^„) -p^(r)ei(p„) -p„(r)ei(p„)]. 



(5) 



M 




(6) 



1=1 
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In turn, the atomic wave functions !f{r — Rj) can be written as a linear superposition of the 
basis functions Xi 

V9(r-R,)=Ec.,X^(r-R.). (7) 

where Rj)} are usually the gaussian basis functions, and the coefficients c^i can be 

obtained from the variational method. 

Then the electronic density of the molecule is, 

p(r) = |^(r)|2 = 5]rf,,^^,x,(r-R.)x,(r-R,), (8) 

ijkl 

where d^kji = 2Cj^Cj;, which can be divided into two parts, namely the net {i = j) and overlap 
7^ j) populations. The latter cannot be ignored in the presence of strong intramolecular 
covalency. Therefore, separating p(r) into its atomic constituents is to split the overlap 
population. However, the way to achieve that is not unique. For example, we can introduce 
a set of weights Wijki due to different criteria such that 

dik,ik = dik,ik + '^ijkidikji / Xfc(r-Ri)x/(r-Rj)rfr, 
^jiji = J2 ('^ - '^ijki)d,kji J Xfc(r-Ri)x;(r-Rj)dr, (9) 

dik,u ~ d^i^^ii, 
then we can rewrite Eq. (jH)) as following 

p(r) ^ ^p,(r) = ^4^^,Xfc(r-R.)x,(r-R,), (10) 

i ikl 

where Pi(r) is the atomic density of atom i. 

In our previous studies, the overlap electronic density is equally separated, i.e., w^ki = 1/2 
in Eq. ©, similar to the MuUiken population analysis (MPA).iii'i2iiLi2i^ In Table HI we 
present the electronic multipole moments of SnClg", ScFg~, Si04~, Se04~, S04^, CIO J, 
C03^, and NO3 calculated from using the MPA, and compare them with the ab initio 
values. We note that for these symmetrical molecules, the MPA preserves total charge and 
zero dipole moment. However, for a molecular ion like V-shaped NO2 or linear CN~, the 
MuUiken population seems inadequate. Given total charge and the dipole moment, 0.26 
a.u., of NO2 obtained from the ab initio calculations (Table H}, a population analysis which 
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TABLE I: Electronic multipole moments of molecule (ABn) calculated from the Mulliken pop- 
ulation analysis. The ab initio values are shown in parentheses. All quantities are in atomic 
units. 







?? 


zz 




^, 


!ZZZ 


CN- 


0.14 (0.17) 


-27.35 


;-29.41) 


7.99 (10.68) 


-143 


(-166) 


NO2- 


0.57 (0.26) 


-19.64 


;-21.64) 


-2.60 (-5.87) 


-60 


(-76) 


NO3- 


0(0) 


-15 


;-i6) 


0(0) 


-38 


(-41) 


cor 


0(0) 


-17 


;-i7) 


0(0) 


-48 


(-48) 


C104 


0(0) 


-111 


;-iii) 


0(0) 


-548 


(-548) 


sof- 


0(0) 


-119 


;-i2o) 


0(0) 


-620 


(-627) 




0(0) 


-143 


;-i43) 


0(0) 


-823 


(-816) 


SiOf- 


0(0) 


-157 


;-i58) 


0(0) 


-1011 


(-1019) 


ScF3- 


0(0) 


-319 


;-3i8) 


0(0) 


-5112 


(-5071) 


SnCir 


0(0) 


-845 


;-844) 


0(0) 


-19834 


(-19668) 



"In the HF calculations, basis set D95* were used for AB and AB2, 6-31G* for AB3 and AB4, and 3-21G* 
for ABq. 

''The electrostatic moments /i (dipole), d (quadrupole), fl (octapole), and <i> (hexadecapole) refer to the 
center of mass of the molecule with the standard orientation defined in GAUSSIAN 98 fRef. l28|) . 

preserves these vales would give rise to -0.092e on N and -0.454e on O. Whereas, using the 
MPA {wijki = 1/2), the charges on the N and O atoms are 0.1624e and -0.5812e, respectively, 
which renders the dipole moment 0.57 atomic unit, overestimated by 120%. 

Therefore, it is desirable to determine Wijki in such a way that the calculated multipole 
moments of the molecule are consistent with the ab initio values. One possible way is to 
evaluate Wijki by fitting the ab initio charge density, as shown in Eqs. © and (fTIH) . with the 
values of multipole moments as constraints. An alternative way is to directly fit the charge 
density, Eq. (fTUj) . with d^i^ji being the parameters and R.j)Xi(r— Rj) being the depen- 

dent variables. To simplify the computation, only the radial parts of x^(r— Rj)x;(r— Rj) 
were kept. This fitting population analysis (FPA) is similar to that proposed by Parker and 
his co-workers as an alternative implementation of the GK model.— 
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C. Free ion and intramolecular interactions 

GAUSSIAN 98 program package^ is employed to obtain the self-consistent solutions to 
the Hartree-Fock-Roothan equations. The atomic orbital basis sets used are a double-zeta 
basis with polarization functions (D95*) for the nitrogen and oxygen atoms. It was shown 
that this choice of basis set could give a good description of free NO^.— As for the sodium 
atoms, we used both the standard 6-31G* basis and the Slater-type orbitals for Na+ taken 
from the Clementi and Roetti table, it turned out the difference between them is small. 

We first consider the free ion case and will discuss the in-crystal ions in Section IV. The 
molecular geometry of NO2 can be described by internal coordinates (-Rnoi; -Rno2 5^ono) 
where -Rnoi and -Rno2 the lengths of the two N-0 bonds, respectively, and 6'ono is 
the angle formed by the two N-0 bonds. The minimum energy geometry of NO^ using 
our basis yields R^q = 1.233 A and 6*01^0 = 116.6°. These structural parameters are 
comparable to the experimental geometry of NO2 in the ferroelectric phase of NaN02,— 
in which R^q = 1.236 A and Oq^^q = 115.4°. The internal vibration frequencies of the 
NO^ group are then calculated at the optimized geometry. Note that the lowest vibration 
frequency (1192 K) of NO2 is considerably higher than the highest libration frequency (318 
K) obtained from Raman spectroscopy^ as well as the order-disorder transition temperature 
(437 K). Therefore, it is justified to treat the internal motion of the nuclei in the NO2 
group within the harmonic approximation, or even as a rigid rotor. To further support this 
argument, we found that only small changes in the geometry of NO^ were observed in the 
high temperature paraelectric phase (-Rno — 1-239 A and ^qno — 114.5°).-^ 

The intramolecular interactions are then represented by the following force field obtained 
from frequency analysis in GAUSSIAN 98:— 

U = Uq + -^no[(-^NOi ~ -^No)^ + (-Rn02 ~ -^No)^] + 2^0No(^ONO ~ ^ONo)^ + 

^2(-^NOi ~ -Rno)(-^N02 ~ -^No) + ^sl-^NOi + -^N02 ~ 2-Rno)(^ONO ~ ^ONo) i^^) 

where Uq = -204.121, R^q = 2.330, O^q^q = 2.036, k^o = 0.728, ^qno = 0-677, = 0.174, 
fcg = 0.066 in atomic units. 

Note that the polarizability of NO2 at its optimized geometry is highly anisotropic, that 
is, with axx = 7.820, ayy = 10.823, azz = 23.825 in atomic units (see Fig. [T] for coordinate 
convention). Thus, one would expect this polarization to seriously affect the intermolecular 
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pair potentials, and thus renders the rigid-ion approximation in question. Nevertheless, after 
extensive ab initio HF calculations with various configurations, which will be presented in the 
next subsection, we find that the intermolecular potentials for the NO2 iNa"*" and NO^:NO^ 
dimers can still be correctly reproduced. 

III. INTERMOLECULAR POTENTIALS 

We scan the potential energy surface of NO^iNa"*" and N0^:N02 obtained from the 
Hartree-Fock method. In these calculations, the geometrical variables of NO^ are frozen at 
their equilibrium values, since we showed previously that the NO2 group in NaN02 could be 
treated as a rigid rotor. However in the ab initio HF calculations, the electronic structure is 
allowed to vary in order to minimize the total energy, thus the electronic polarization effects 
are included. 

In our calculations of intermolecular interactions, one NO2 is fixed with its center of mass 
being the origin of the coordinate system, the dipole vector pointing to the y axis and the 
0-0 line being aligned parallel to the z axis (see Fig.Q). Then, Cartesian coordinates (x, y, 
z) are the position of Na"*" or the center of mass of another NO2 . To describe the orientation 
of the NO2" molecule at (x, y, z), we use the angles a and [3 of the dipole moment vector of 
the NO^ molecule and an angle 7 involving the rotation of NO^ around its dipole vector. 

In order to study the effects of the different partition schemes, mentioned in Section IIB, 
on the pair potentials, we performed three different calculations, namely (i) MPA with pair 
potentials [Eq. (0)], (ii) FPA with pair potentials, (iii) FPA with non-pair potentials as 
shown in Eq. We shall refer them as models I, II, and III, respectively, from now on. 
Recall that model I overestimates the dipole moment of NO^ as mentioned in Section IIB. 
Furthermore, we shall show in the following that the electronic polarization effect could be 
revealed from examining the differences between model I and II, while the errors due to the 
pairwise additive approximation employed by the GK model could be analyzed from the 
differences between model II and III. 

In Figs. 121 and El we compare the pair potentials obtained from these models and those 
from the ab initio HF calculations. In Fig. |2lwe show the results for the N02:Na"*' dimer: 
we find that in both models I and II, the overall shapes of the GK potentials as a function of 
molecular separation are in agreement with the ab initio results, particularly in Fig.|2fd). On 
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FIG. 2: NO^-Na"*" intermolecular potential energy curves as a function of R for various configura- 
tions: (0,i?, 0), 0,0), (0, — i?, 0), and (0,0, i?), where {x,y,z) is the location of Na+. Different 
lines represent the ab initio HF model (solid), model I (dashed), model II (dotted), and model III 
(dashed and dotted). 

the other hand, the electronic polarization effect also manifests itself in Fig. |2] based on the 
following two observations. First, notice Figs.Efb) and|2fd), models I and II fit the ab initio 
results through the entire range in {R, 0, 0), whereas in (0, 0, R) models I and II overestimate 
the potentials when 5 < i? < 8. We attribute that to the anisotropic polarizability of NO2 
{a^x < Ciyy < Cizz)i thus the electronic cloud of NO2 is most unlikely to be polarized along 
the (-R, 0, 0) direction, along which the polarization is the weakest, rather than the (0, 0, R) 
direction. Second, in Fig. |2fa), (0, i?, 0), the minimum potential energy in model II is closer 
to the ah initio values than model I, whereas in (0, — i?, 0), where Na^ is located below 
NO^, we found the situation reversed. To understand this, we observe that in the (0, R, 0) 
configuration, Na"*" is closer to N than O, thus in the ah initio calculations the electrons were 
attracted toward the N atom and thus led to a smaller dipole moment. Therefore, model II, 
which produced smaller dipole moment than model I, tends to get closer to the HF results. 
Obviously, the results will be opposite if the situation is reversed. 

Among the four different configurations. Figs. Ufa) through |2td), and methods except 
model III, the lowest N02-Na^ potential energy takes place at (0, — i?, 0), which will be 
shown later to give rise to the lowest energy structure of NaN02. Both models II and III use 
FPA, the only difference between them is that model II employs the pair potentials. It seems 
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FIG. 3: NO2 -NO2 intermolecular potential energy curves as a function of rotation angle. The 
N02^ molecules are initially in parallel alignment at separation (1.82 A, 2.83 A, 2.69 A) and then 
one of them rotates around one of the (a) x, (b) y and (c) z axes through its center of mass. 

that in model II the errors caused by the pairwise additive approximation is compensated 
by the errors due to FPA. 

Similarly, we show in Fig. IHlthe N02-N0^ intermolecular potentials using the same four 
methods. The configurations are chosen as follows: the two NO2 molecules are initially 
parallel at their experimental low-temperature separation (1.82 A, 2.83 A, 2.69 A) and then 
one of them rotates around each of the x, y, z axes through its center of mass, as shown in 
Figs. E^a) through infc). The results of models I-III agree reasonably well with the ah initio 
calculations. On closer examination. Models II and III fit better to the ah initio results 
than I. We believe this is due to the fact that the two nitrite ions are separated so far that 
the main contribution to the intermolecular potential is mostly electrostatic. Based on the 
experience acquired in Section IIB, FPA usually gives a better description of the rigid NO2 
molecule than MPA, hence the results obtained from models II and III which employ the 
FPA are considerably better. In addition, the short range interaction is also small, thus the 
difference between model II and III is rather small. 

Summarizing this section, in spite of the presence of electronic polarization when two 
molecules are brought closer, the intermolecular potentials for the NO2 :Na+ and N0^:N02 
dimers could be correctly reproduced using models I and II. Furthermore, model II appears 
to give better description of NO2 than model I as shown in Fig. El particularly Figs. Efa) 
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TABLE II: Experimental and theoretical structural parameters for the Im2m phase (III) of NaN02. 
Lattice constants are given in A. 



Parameters 


Experiment" 


Model I 


Model II 


Model III 


a 


3.5180 


3.3889 


3.5013 


3.7353 


b 


5.5350 


5.4542 


5.5485 


5.4257 


c 


5.3820 


4.9254 


4.8403 


4.9669 


y/b of N (2a) 


0.5437 


0.0498 


0.0433 


0.0586 


y/b of Na (2a) 


0.0781 


0.5537 


0.5492 


0.5437 


y/b of (4d) 


-0.0443 


-0.0704 


-0.0740 


-0.0610 


z/c of (4d) 


0.1962 


0.2111 


0.2151 


0.2098 



Trom X-ray diffraction experiments at 120 K, see Ref. 
and Of c). 
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IV. MOLECULAR DYNAMICS SIMULATIONS 



After the potential energy surface for NaN02 has been obtained, we are prepared to un- 
dertake MD simulations. Long-range Coulombic interaction in the crystal is represented by 
electrostatic interaction among point charges calculated from the population analysis, while 
each of the short-range pair potentials are fitted by using an exponential-polynomial func- 
tion accurate within 0.1%.— In the following description, the x, ?/, z directions correspond 
to the crystallographic a, 6, c directions of NaN02, respectively, see Fig. [TJ 



A. Lattice relaxation 



Before we proceed with the molecular dynamics simulations, we perform lattice relaxation 
for the ferroelectric structure of NaN02 both with and without the Ivrtlm space group 
symmetry constraints. This relaxation procedure provides the crystal structure with zero 
force on each atom, that is an energy extremum; it also produces a test to the PES because 
the resultant structures have to agree reasonably with the experimental data for further 
simulations to be reliable. We perform both static and dynamic relaxations: the static 
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one is an application of the Newton-Raphson algorithm and usually results in finding a 
local minimum of the energy, and the dynamic one is a simulated annealing calculation for 
overcoming that limitation. We start the static lattice relaxation with the experimental 
parameters. In Table [H] we present the lattice and basis parameters deduced from the 
experiments and static relaxation. In all cases, the static relaxation produced essentially 
the same structure that strongly resembles the experimental structure. Most of the lattice 
constants in the relaxed structure are shorter than the experimental values (by 3.7%, 1.5%, 
and 8.5% for a, b, and c, respectively, in model I, and by 0.5%, -2.4%, and 10% for a, b, and 
c, respectively, in model II). Hence the calculated volume is smaller than the experimental 
one by 13% for model I and 10% for model II, a common feature for simulations using the 
GK model, which will be addressed in more detail in the next subsection. 

Next, we go on to relax the statically relaxed crystal structure to zero temperature using a 
simulated annealing algorithm, in which the amount of kinetic energy in the molecules slowly 
decreases over the course of the simulation. We find that the (zero temperature) ground 
states in models I and II are close to the statically relaxed structures, whereas there are 
substantial changes taking place in model III. By monitoring the orientations of the nitrite 
ions, we find that the ground structure in model III, still orthorhombic with a = 3.90494 A, 
b = 4.8441 A, and c = 5, 0770 A, is ferroelectric with the dipole moments of NO^ aligned 
along the a axis rather than the experimental b axis. So we conclude that the PES given by 
models III did not refiect reality. This concurs with the previous discussion (Section III) on 
the intermolecular potentials. In the following we use only models I and II to simulate the 
phase transition in NaN02. 

B. MD simulations 

Using the isoenthalpic, isobaric ensemble, our MD simulation is started with a zero- 
temperature zero-pressure orthorhombic cell (4a x 46 x 4c), consisting of 512 atoms. Periodic 
boundary conditions are imposed to simulate an infinite crystal. The MD calculations are 
carried out in the Parrinello- Rahman scheme^ which allows both the volume and the shape 
of the MD cell to vary with time. The calculation of the energies and forces was handled by 
the Ewald method. A time step of 0.002 ps was used to integrate the equations of motion. 
In our heating runs, we raise the temperature of the sample in stages, 20 K each time, up 
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FIG. 4: Temperature variation of lattice constants o, b, c (solid, dashed, dotted lines, respectively; 
left scale) and volume of the unit cell (open circles; right scale) for (a) model I and (b) model II. 

to 1000 K. At each stage, the first 2000 time steps were employed to equilibrate the system, 
then 10000 time steps were collected for subsequent statistical analysis. Since our simulations 
employ periodic boundary conditions, we cannot distinguish the incommensurate structure 
(i.e., phase II of solid NaN02). 

In Figs.|l]through[3 we demonstrate that as the MD cell is heated, it undergoes two phase 
transitions: in the first one, the system retains its orthorhombic structure with a change 
of space group from Im2m to Immm, in agreement with the experiments. The critical 
temperature Tq is around 400 K for model I and 313 K for model II, respectively. In the 
second transition, the crystal structure violently changes from orthorhombic to tetragonal 
at temperature (Tm) which is around 500 K for model I and 350 K for model II, as shown in 
Fig. 13 However, we argue that the crystal has already melted before this type of transition 
could be observed in reality. 

To investigate the mechanism of the polarization reversal of NO^, we monitor the crystal 
polarization and display the results in Fig.^l Let the dipole moment of anion i be nij and the 
quadrupole moment be t?, calculated by using the point charges on the and O atoms. Then 
the mean dipole moment per anion at temperature T is M(T) = J2i {^i) where N = 128 
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FIG. 5: Mean dipole moment Mh{T) and quadrupole moment Q of the whole NaN02 crystal as a 
function of temperature for the MD runs for (a) model I and (b) model II. 

is the number of NO2 in the MD cell and the brackets denote an average over the MD run. 
In addition, we define the antiferroelectric polarization as Q = X^i Gxp(— tt ■ Rj) (nij) /N 
where tt = (tt, tt, tt) and Rj is the lattice vectors associated with the i-th ion. Within 
our statistical uncertainty we find over all temperature range = = Q = 0, while 
My{T < Tq) > and My{T > Tq) = 0. This fact confirms that the transition taking place 
at Tq is the ferroelectric-paraelectric phase transition. Furthermore, we calculated the mean 
quadrupole moment = J2i (^j) When the dipole vector of a NO2 is aligned along the 
b axis, ~ 0.00, ~ -0.04, ~ -4.49 for model I and ~ 0.00, ^yy ~ -0.20, 
— —3.52 for model II; thus {Qxx + Qyy)/'2'dzz ^ 1- This relation holds as the NO2 
ion rotates around the c axis; nevertheless, one would expect {Qxx + Qyy)/'^Qzz = 1 when 
the NO^ ion rotates without directional preference. The fact that {Qxx + Qyy)/'^Qzz <^ 1 
for T < Tm (Fig. E]) reveals that the NO^ anions rotate primarily about the c axis. When 
T > Tjn, {Qxx + ^yy)/'^^zz — 1, i-©., NaN02 becomes an orientational liquid. 
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FIG. 6: Diagonal elements of the mean-square atomic displacements Ua vs temperature, (a) Na, 
(b) N, and (c) O atom. 

Further, in Fig. IHl we show the mean-square atomic displacements Uu = {uj) where 
i = 1, 2, 3 denotes the displacements along the a, b, c axes, respectively. Different models 
of NO^ reversal are expected to satisfy the following conditions: (1) rotation around the 
c axis: t/22(N), f/33(N) < Uii{N) and U22{0) , UssiO) < f/ii(0); (2) rotation around the 
a axis: ?7ii(N), f/22(N) < U^siN) and f/ii(0), f/33(0) < f/22(0). This figure relates to 
recent X-ray experiments which used the same quantities to investigate the polarization 
reversal mechanism.— The experiments confirmed that the first condition holds for both 
ferroelectric and paraelectric phases. Another important feature revealed by the experiments 
is that f/22(Na), f/33(Na) < f/ii(Na) in the ferroelectric phase, whereas t/ii(Na), t/33(Na) < 
?722(Na) in the paraelectric phase. That is, f/ii(Na) and [/22(Na) are reversed across Tq- 
These features are reproduced in Fig. IHlwith exception of f/ii(0), f/33(0) < f/22(C)) in the 
paraelectric phase. This means the NO2 motions in our simulations are more mobile than 
those in the real crystal, rendering the simulated transition temperatures lower than the 
experimental values of Tc — 437 K and the melting temperature 550 K. In other words, the 
barriers to NO2 rotation in our models are too small. 
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FIG. 7: Atomic positions of NaN02 viewed from the a direction obtained from the MD simulation 
for model I at (a) T = 198 K, (b) T = 320 K, (c) T = 449 K, and (d) T = 535 K. 

In addition, in Fig. [71 we show the average crystal structures of NaN02 at different 
temperatures. The ellipsoids in these pictures represent the root-mean-square deviations 
of the atoms from their average positions and thus indicate the thermal motions of these 
atoms. The c-axis rotation mode can be clearly seen, particularly in Fig. El^c). 

According to the calculations described in the previous section, model II generally gives 
a better description of free NO2 and the intermolecular potential energies than model I. 
However, the simulation based on model I matches closer the experimental results than 
that based on model II, that is, Tq and Tm predicted by model I are closer to experiment, 
this indicates that the crystal fields and polarization effects are particularly important for 
quantitatively studying the NaN02 system, where the ferroelectric structure is considerably 
stabilized by anion polarization effects. Actually, the MPA employed by model I does not 
preserve the ab initio dipole moment of free NO2 . Rather, it overestimates the dipole 
moment by 120%, thus leading to higher NO^ rotational barriers than those predicted by 
model II, which in turn raises the transition temperature and provides a better simulation 
in comparison with the experiments. 

It is worth mentioning the less desirable agreement between theoretical and experimental 
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volumes, namely, the 13% discrepancy for model I and 10% for model II. To address this 
we make one simple change: by following Waldman and Gordon,— we increase the kinetic 
energy term in the Gordon-Kim potentials by 5%, this reduces the discrepancy to 9% for 
model I and 6% for model II. Having done this we rerun the MD to obtain values of of 
360 K for model I and 303 K for model II. While this change worsens the value for model I, 
the value for model II is virtually unchanged. And in both cases the transition mechanism 
is unaltered. Thus the slight hardening of the short-range potentials removes most of the 
volume discrepancies. However, there is no material change in the mechanism of the phase 
transition. This robustness of the results with respect to minor variations in the potential 
demonstrates that our basic conclusion remain valid. 

C. Rotational barriers 

Based on the previous simulation results, the order-disorder phase transition in NaN02 
involves the rotation of the nitrite ions. We devise a scheme to calculate the three rotational 
barriers for NOg around the crystallographic a, b, and c axes with its center of mass fixed: 
we start from the experimental ferroelectric structure^^ as the ground state and calculate 
its energy difference with one of the two nitrite ions in the unit cell being rotated about 
the respective axis. The results are shown in Figs. IHa) and IHl^b) . The bottom of each 
barrier, zero angle, is in the ferroelectric structure. For both models I and II, the rotation 
around the c-axis has an energy barrier 5-10 times smaller than those of the other rotations, 
which is a characteristic of nitrites^^^ Hence, our calculations unambiguously reveal that 
the reorientation of NO^ in the paraelectric phase occur essentially by rotations around the 
c axis. In addition, the barriers calculated in model I are higher than that in model II, 
confirming that the transition temperature in model I will be higher than in model II. 

In Figs. IHl^c) andlHl^d), we also plot the contribution to the rotational barriers purely 
from the electrostatic interaction, that is, in the point-charge model with the short range 
forces deleted. Comparing Fig. |Hl^a) with Fig. |Hfc), or Fig. |Hl^b) with Fig. IHl^d), we notice 
that the point-charge model gives rise to a quite different rotation barrier about the a axis: 
it bottoms at about ±90° and is lower than the ground state energy due to the omission 
of short-range interactions, which comes from overlap of the charge cloud of an atom with 
those of its neighbors. 
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FIG. 8: Rotational barriers of one of the two nitrite ions in the unit ceh around the a, b, and c 
axes with its center of mass fixed. (a)(b) for the GK model, (c)(d) for the point charge model. 
Left and right panels are for model I and II, respectively. 

V. CRYSTAL-FIELD EFFECTS 

In this section, we investigate the crystal-field effects on the NO2 ion, which include 
electrostatic interaction from the background lattice, overlap compression of the NO2 charge 
cloud through interaction with its neighbors, and charge-transfer between molecules which 
is usually small in ionic crystals. In the studies of monatomic iona^^ and cyanides^. Fowler 
et al. showed that these effects could be successfully described by embedding the ion of 
interest, or a cluster consisting of the ion and its first shell of neighbors, into a lattice of 
point charges. 

We therefore perform HF calculations based on the following two models.— First, the 
crystal field of ferroelectric NaN02 is simulated by placing the nitrite ion at the center of a 
4x4x4 orthorhombic point charge lattice with spacings equal to the experimental lattice 
parameters. Charges in the faces of the lattice are scaled to maintain overall neutrality. All 
anions except the central NO2 are represented by single point charges. Hence, there are 
174 point charges surrounding the NO2 ion. Calculations of this type are referred to as 
CRYST. Obviously, in the CRYST calculation we take into account only the crystal-field 
effect arising purely from electrostatic interaction. 

At the next level of sophistication, we replace the six nearest positive charges of the 
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FIG. 9: Dipole moment of the central NO^ in a 4 x 4 x 4 lattice as it rotates around the a, b, c 
axes through the center of its (Na"'')6 cage, (a) CLUST and (b) CRYST. 

central NO2 in the above lattice by the Na cations. Calculations of this type are referred to 
as CLUST. In both CRYST and CLUST calculations, the geometrical structure of the NO2 
is fixed at its experimental values. We employ the same basis set, D95*, for the in-crystal 
NO^ ion as for the free NO^ ion. In order to keep the CLUST calculations to a manageable 
size, we use the minimal basis set, ST0-3G, for the Na^ ions unless specified. The cations, 
however, are relatively insensitive to the crystal environment and they are included here only 
to account for their compressing effect on the NO2 wave functions. We find that adding 
extra basis functions to Na"*" will not change the results significantly. 

The in-crystal NO2 initially points in the b direction as in the ferroelectric phase of 
NaNO#. Its dipole moment is 0.636 (CRYST) and 0.661 (CLUST) Debye, close to that in 
the free ion model, 0.661 Debye. Thus it appears that the crystal-field effect is small for this 
configuration. Subsequently, we rotate the NO2 about the a, b, and c axes and calculate the 
dipole moment of the rotated NO^. The rotation center is taken to be the center of a (Na"'")6 
cage formed by the 6 neighboring sodium ions of the central NO2 (0,0.279 A, 0) in the 
coordinate convention of Fig. ^ The CLUST and CRYST results are depicted in Figs. IHl^a) 
and El^b), respectively. Clearly, the dipole moment of the NO2 changes considerably as 
it rotates, indicating strong crystal field effects on the reorientation of the NO^. Since the 
electron density of the NO2 is compressed by its 6 neighboring sodium cations, the variation 
in magnitude of its dipole moment is smaller in the CLUST model than that in the CRYST 
model. The electron cloud of the NO2 is most and least variable when it rotates around the 
a and b axes, respectively; for the c-axis rotation, which has the lowest rotational barrier. 
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FIG. 10: ab initio barriers to rotation of the central NO^ in a 4 x 4 x 4 lattice around the a, b, c 
axes through the center of its (Na+)6 cage, (a) CLUST and (b) CRYST. 

the dipole moment goes down from 0.661 Debye at 0° to 0.534 Debye at 180° in the CLUST 
model, as shown Fig. IHl^a). 

In the context of population analysis, increase of the dipole moment of NO^ implies 
that more electrons are distributed on the O atom, i.e., electrons are flowing from the 
nitrogen atom to the oxygen atoms. Conversely, decrease of the dipole moment indicates a 
reversal in electron transfer. Therefore, we have demonstrated considerable intramolecular 
charge-transfer, although the intermolecular charge-transfer is usually small in ionic crystals. 
This intramolecular charge-transfer could be further elucidated based on the language of 
molecular orbitals (MOs): in the minimal basis set of atomic orbitals (AOs) on nitrogen 
and oxygen, each atom of the NO2 molecule contributes one p-orbital perpendicular to the 
molecular plane. Thus their linear combinations, which are determined by the crystal fleld, 
form three different vr MOs extended over the entire molecule, thus leading to the above- 
mentioned intramolecular charge-transfer. Note that the mean HF polarizability of free 
NO^, 14.156 in the atomic units, is much higher than that of free Na^, 0.343, calculated 
by using the 6-31* basis. Therefore, it is reasonable to expect that NO2 in solid NaN02 
encountering different crystal-field environments as it rotates, redistributes its charge among 
the three constituent atoms to lower its energy. 

In Fig. ^Jwe show that the rotational barriers of the central N02^. On the whole. The 
CLUST results are similar to those in Figs. Efa) and|Hl^b), while the CRYST results are 
similar to those in Figs. |Hfc) andlH^d). The reason is that in the CRYST calculations we 
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FIG. 11: In the CLUST model, NO^ rotates around the a, b, c axes through its center of mass, (a) 
dipole moment of NO2 , and (b) rotational barriers. 

consider only the crystal-field effects originating purely from electrostatic interactions with 
the background point charges, similar to the point charge model used to obtain Figs. |H^c) 
and Efd). Whereas, in the CLUST calculations, overlap compression is also taken into 
account, which is reflected in Figs. Efa) and|Hfb) as inclusion of short range repulsion. In 
Fig. lior a). the barrier to the c-axis rotation is the lowest, but comparable with that to the b- 
axis rotation. This feature is caused by the fact that all background anions are represented 
by single point charges in our CLUST and CRYST calculations; however, we anticipate 
that restoring multipole moments of these background anions would increase the barrier 
difference among rotations about the a, b, and c axes. 

To further elaborate the electronic polarization effect on N02^ arising from the crystal 
environment, we change the rotation center to the center of mass of the NO2 which was 
assumed in the model by Kremer and Siems.— In this case, the dipole moment of NO^ 
and the rotational barriers are presented in Figs. [Ufa) and lTTT b). respectively. The dipole 
moment changes in a different way from that shown in Fig. [Ufa). In particular, the dipole 
moment in the ferroelectric phase (0°-rotation) is larger than at 180°-rotation in Fig. El^a), 
while it is smaller in Fig. ITlT a). On the other hand, there is no qualitative discrepancy 
between Fig. lTlT b) and Fig. ITUT a): the main differences are: the rotational barrier about the 
a axis has risen by 135% while the barrier about the c axis is depressed by 26%. This means 
that the order of rotational barriers is enhanced by the change of the rotation center. 

Although strong crystal field effects have been revealed by these ab initio calculations, 
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the rotational barriers obtained from the polarizable-ion models are in qualitative agreement 
with those from the rigid-ion models, confirming that the rigid-ion model is capable of 
describing the phase behavior in NaN02. 

VI. CONCLUDING REMARKS 

We have presented MD simulations of NaN02 using a hybrid a priori method consist- 
ing of ah initio calculations and Gordon-Kim electron gas theory to analytically calculate 
the crystal potential surface. This method has been carefully examined by using differ- 
ent population analysis methods. We have carried out ah initio Hartree-Fock calculations 
of the intermolecular interactions for NO2 :Na^ and NO2 :N0^ dimers and concluded that 
the pair potentials of the rigid-ion model can correctly reproduce the ah initio results. We 
demonstrated that a rigid-ion model is capable of describing phase behavior in solid NaN02. 

The crystal-field effects on the NO2 ion are also addressed in two polarizable-ion models 
for which the ferroelectric phase of NaN02 was found to have a larger dipole moment of 
NO^ than the 180°-rotation phase. Remarkable internal charge-transfer effect is found 
to be stabilizing the crystal structure of NaN02. In our MD simulations, two rigid- ion 
models using MPA and FPA, respectively, have been studied. The one using MPA, which 
overestimates the dipole moment of NO2 , gives rise to the more comparable results with the 
experiments, since such overestimation also stabilizes the crystal structure, thus mimics the 
anion polarization effect. To quantitatively simulate NaN02, a more elaborate polarizable- 
ion model is needed. 
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